7.2.1 生成坡向坡度数据

## 基于投影坐标系绘制坡向和坡度：已经进行二维地图优化高程显示：
higground = dem >600
slope = terrain(dem,opt ='slope')
aspect  = terrain(dem,opt='aspect')
hill = hillShade(slope,aspect)
plot(hill,col= grey(0:100/100),legend =F)
plot(dem,col =ranibow(25,alpha=0.35,add=T)

7.2.2 rayshader构建三维立体地图

## 构建三维地图：
## 凡是三维的数据都可以考虑使用这个R包，能做一些相对有趣的图和可视化工作；
# 最终可以实现的效果：3点可视化物种分布点及物种分布图；

library(rayshader)
setwd("E:/环境数据/dem/")
#Here, I load a map with the raster package.
loadzip = tempfile()
download.file("https://tylermw.com/data/dem_01.tif.zip", loadzip)
localtif = raster::raster( "./dem_01.tif")


plot(localtif)

#And convert it to a matrix:
elmat = raster_to_matrix(localtif)

#We use another one of rayshader's built-in textures:
##  构建2d可视化：
elmat %>%
  sphere_shade(texture = "desert") %>%
  plot_map()
# 更改太阳倾角：
elmat %>%
  sphere_shade(sunangle = 45, texture = "desert") %>%
  plot_map()

# 基于高程构建水流环境：
elmat %>%
  sphere_shade(texture = "desert") %>%
  add_water(detect_water(elmat), color = "desert") %>%
  plot_map()

## 颜色变暗：
elmat %>%
  sphere_shade(texture = "desert") %>%
  add_water(detect_water(elmat), color = "desert") %>%
  add_shadow(ray_shade(elmat), 0.5) %>%
  plot_map()

## 添加大气散射光环境层：
elmat %>%
  sphere_shade(texture = "desert") %>%
  add_water(detect_water(elmat), color = "desert") %>%
  add_shadow(ray_shade(elmat), 0.5) %>%
  add_shadow(ambient_shade(elmat), 0) %>%
  plot_map()

## plot_3d中zscale的值增大会导致山峰变得更加锐利：默认是10
elmat %>%
  sphere_shade(texture = "desert") %>%
  add_water(detect_water(elmat), color = "desert") %>%
  add_shadow(ray_shade(elmat, zscale = 3), 0.5) %>%
  add_shadow(ambient_shade(elmat), 0) %>%
  plot_3d(elmat, zscale = 10, fov = 0, theta = 135, zoom = 0.75, phi = 45, windowsize = c(1000, 800))
# 使用下面这个代码保存截图：
Sys.sleep(0.2)
render_snapshot()


## 添加比例尺和指南针：
elmat %>%
  sphere_shade(texture = "desert") %>%
  add_water(detect_water(elmat), color = "desert") %>%
  add_shadow(ray_shade(elmat, zscale = 3), 0.5) %>%
  add_shadow(ambient_shade(elmat), 0) %>%
  plot_3d(elmat, zscale = 10, fov = 0, theta = 135, zoom = 0.75, phi = 45, windowsize = c(1000, 800))

render_camera(fov = 0, theta = 60, zoom = 0.75, phi = 45)
render_scalebar(limits=c(0, 5, 10),label_unit = "km",position = "W", y=50,
                scale_length = c(0.33,1))
render_compass(position = "E")
render_snapshot(clear=TRUE)

## 深度渲染可视化：耗费时间的呀！而且不怎么好看！
elmat %>%
  sphere_shade(texture = "desert") %>%
  add_water(detect_water(elmat), color = "desert") %>%
  plot_3d(elmat, zscale = 10, fov = 0, theta = 60, zoom = 0.75, phi = 45, windowsize = c(1000, 800))

render_scalebar(limits=c(0, 5, 10),label_unit = "km",position = "W", y=50,
                scale_length = c(0.33,1))

render_compass(position = "E")
Sys.sleep(0.2)
library(rayrender)
render_highquality(samples=200, scale_text_size = 24,clear=TRUE)

## 可视化水层：并添加水的颜色：
montshadow = ray_shade(montereybay, zscale = 50, lambert = FALSE)
montamb = ambient_shade(montereybay, zscale = 50)
montereybay %>%
  sphere_shade(zscale = 10, texture = "imhof1") %>%
  add_shadow(montshadow, 0.5) %>%
  add_shadow(montamb, 0) %>%
  plot_3d(montereybay, zscale = 50, fov = 0, theta = -45, phi = 45,
          windowsize = c(500, 300), zoom = 0.75,
          water = TRUE, waterdepth = 0, wateralpha = 0.5, watercolor = "lightblue",
          waterlinecolor = "white", waterlinealpha = 0.5)
Sys.sleep(0.2)
render_snapshot(clear=TRUE)


##  baseshape = "circle"修改成图的形状；circle为圆形，hex为五角形；
montereybay %>%
  sphere_shade(zscale = 10, texture = "imhof1") %>%
  add_shadow(montshadow, 0.5) %>%
  add_shadow(montamb, 0) %>%
  plot_3d(montereybay, zscale = 50, fov = 0, theta = -45, phi = 45, windowsize = c(1000, 800), zoom = 0.6,
          water = TRUE, waterdepth = 0, wateralpha = 0.5, watercolor = "lightblue",
          waterlinecolor = "white", waterlinealpha = 0.5, baseshape = "circle")

render_snapshot(clear = TRUE)


## 使用render_label()功能可以添加文本标签，该功能还允许您自定义线条类型，颜色和大小以及字体：
montereybay %>%
  sphere_shade(zscale = 10, texture = "imhof1") %>%
  add_shadow(montshadow, 0.5) %>%
  add_shadow(montamb,0) %>%
  plot_3d(montereybay, zscale = 50, fov = 0, theta = -100, phi = 30, windowsize = c(1000, 800), zoom = 0.6,
          water = TRUE, waterdepth = 0, waterlinecolor = "white", waterlinealpha = 0.5,
          wateralpha = 0.5, watercolor = "lightblue")
render_label(montereybay, x = 350, y = 160, z = 1000, zscale = 50,
             text = "Moss Landing", textsize = 2, linewidth = 5)
render_label(montereybay, x = 220, y = 70, z = 7000, zscale = 50,
             text = "Santa Cruz", textcolor = "darkred", linecolor = "darkred",
             textsize = 2, linewidth = 5)
render_label(montereybay, x = 300, y = 270, z = 4000, zscale = 50,
             text = "Monterey", dashed = TRUE, textsize = 2, linewidth = 5)
render_label(montereybay, x = 50, y = 270, z = 1000, zscale = 50,  textcolor = "white", linecolor = "white",
             text = "Monterey Canyon", relativez = FALSE, textsize = 2, linewidth = 5)


## 添加路径、点和shp，使用sf包；
## 添加polygons：
library(sf)
montereybay %>%
  sphere_shade(texture = "desert") %>%
  add_shadow(ray_shade(montereybay,zscale=50)) %>%
  plot_3d(montereybay,water=TRUE, windowsize=c(1000,800), watercolor="dodgerblue")
render_camera(theta=-60,  phi=60, zoom = 0.85, fov=30)

#We will apply a negative buffer to create space between adjacent polygons:

mont_county_buff = sf::st_simplify(sf::st_buffer(monterey_counties_sf,-0.003), dTolerance=0.001)

render_polygons(mont_county_buff,
                extent = attr(montereybay,"extent"), data_column_top = "ALAND",
                scale_data = 300/(2.6E9), color="chartreuse4",
                parallel=TRUE)
render_highquality(clamp_value=10,sample_method="stratified")


## 添加点：
montereybay %>%
  sphere_shade(texture = "desert") %>%
  add_shadow(ray_shade(montereybay,zscale=50)) %>%
  plot_3d(montereybay,water=TRUE, windowsize=c(1000,800), watercolor="dodgerblue")
render_camera(theta=-60,  phi=60, zoom = 0.85, fov=30)

moss_landing_coord = c(36.806807, -121.793332)
x_vel_out = -0.001 + rnorm(1000)[1:500]/1000
y_vel_out = rnorm(1000)[1:500]/200
z_out = c(seq(0,2000,length.out = 180), seq(2000,0,length.out=10),
          seq(0,2000,length.out = 100), seq(2000,0,length.out=10))

bird_track_lat = list()
bird_track_long = list()
bird_track_lat[[1]] = moss_landing_coord[1]
bird_track_long[[1]] = moss_landing_coord[2]

for(i in 2:500) {
  bird_track_lat[[i]] = bird_track_lat[[i-1]] + y_vel_out[i]
  bird_track_long[[i]] = bird_track_long[[i-1]] + x_vel_out[i]
}

render_points(extent = attr(montereybay,"extent"),
              lat = unlist(bird_track_lat), long = unlist(bird_track_long),
              altitude = z_out, zscale=50, color="red")
render_highquality(point_radius = 1,sample_method="stratified")



#Combining base R plotting with rayshader's spherical color mapping and raytracing:
## Not run: 
##  绘制地形阴影：
montereybay %>%
  sphere_shade() %>%
  add_overlay(height_shade(montereybay),alphalayer = 0.6)  %>%
  add_shadow(ray_shade(montereybay,zscale=50)) %>%
  plot_map()

## End(Not run)

#Add contours with `generate_contour_overlay()`
# 添加等高线：
montereybay %>%
  height_shade() %>%
  add_overlay(generate_contour_overlay(montereybay))  %>%
  add_shadow(ray_shade(montereybay,zscale=50)) %>%
  # plot_map()
  plot_3d(montereybay, zscale = 30, fov = 0, theta = 135, zoom = 0.75, phi = 45, windowsize = c(1000, 800))



### 使用自己的数据进行模拟实验：####
## 注意这里构建dem转投影坐标系时，需要使用双插值方法进行；
setwd("E:/环境数据/dem/")
dem <- raster("./wc2.1_30s_elev.tif")
source("E:/sjdata/code/sdmenm.R")
all <- read.csv("E:/sjdata/F1_SP/Hgya3.csv")[,2:3]
range <- mcp(all) %>% rgeos::gBuffer(width =1)
dems <- mask(crop(dem,range),range)
## 将地理坐标系转为投影坐标系：

demms<- projectRaster(dems, crs = CRS("+init=epsg:3857") )
plot(demms)
## 数据转为矩阵：
demmss = raster_to_matrix(demms)
## 重要的一步：将矩阵中的na值等于0
demmss[is.na(demmss)] <- 0


## 添加点：
demmss %>%
  sphere_shade(texture = "desert") %>%
  add_water(detect_water(demmss, min_area = 15,zscale=10), color = "desert") %>%
  # add_shadow(ray_shade(demmss, zscale = 3), 0.5) %>%
  add_shadow(ambient_shade(demmss), 0) %>%
  ## z值低于50会过于尖锐；
  plot_3d(demmss, zscale = 80, fov = 0, theta = 135, zoom = 0.75, phi = 45, windowsize = c(500, 300))

render_points(extent = attr(dems ,"extent"),
              lat = unlist(all$latitude),long =  unlist(all$longitude),
              altitude = z_out, zscale=20, color="red")
render_water(demmss,zscale=10)

render_camera(fov = 0, theta = 60, zoom = 0.75, phi = 45)
render_scalebar(limits=c(0, 250, 500),label_unit = "km",position = "W", y=50,
                scale_length = 1)
render_compass(position = "N",scale_distance =1.1, compass_radius=80)
render_snapshot()

extent(demms)
(demms@extent@ymax - demms@extent@ymin ) /1000

## 尝试在海拔的基础上添加温度图层：####
bio1 <- raster("E:/环境数据/wc2.0_30s_bio/wc2.0_bio_30s_01.tif")
bio <- mask(crop(bio1,dems),dems)

plot(bio)
biox<- projectRaster(bio, crs = CRS("+init=epsg:3857") )
plot(demms)
## 数据转为矩阵：
bioxc = raster_to_matrix(biox)
## 重要的一步：将矩阵中的na值等于0
bioxc[is.na(bioxc)] <- 0


## 一种可行的实现方法：
landsat_mat_list <- lapply(as.list(landsat), function(x){t(raster_to_matrix(x)/255)})
landsat_mat_list[[4]] <- t(raster_to_matrix(raster(nrows = nrow(landsat), ncols = ncol(landsat), ext = extent(landsat), resolution = raster::res(landsat), vals = 0.9)))

library(abind)
landsat_mat <- do.call(abind,list(landsat_mat_list,along = 3))




hill <- demms %>%
  sphere_shade() 

hill <- add_overlay(hill, landsat_mat)

rgl::clear3d()

plot_3d(hillshade = hill, heightmap = dem_mat, 
        windowsize = c(1000, 600), zscale = zscale, 
        theta = 160, zoom = 0.5, phi = 35, baseshape = "circle")



## 不能使用library(rayshader)这个包，需要更换：
library(rasterVis)
## 添加投影坐标系下的栅格；前置高程，后置叠加栅格；
plot3D(demms,drape= biox , adjust = FALSE)

7.2.3 plot3d

spatialEco::tri()也同
## 样可以计算；
## 快速地形的3d可视化：
library(rasterVis)
plot3D(valley)  ## 其中vally是对应的dem
这个可以使用R包的raster进行计算（terrain()）；
slope <- terrain(dem, "slope")
aspect <- terrain(dem, "aspect")
TPI <- terrain(dem, "TPI") # Topographic Position Index
TRI <- terrain(dem, "TRI") # Terrain Ruggedness Index
roughness <- terrain(dem, "roughness")
flowdir <- terrain(dem, "flowdir")